The pyActigraphy package implements several rest/activity detection algorithms.
Broadly speaking, they can be divided in two categories:
Epoch-by-epoch rest/activity scoring (inspired by PSG):
Cole-Kripke’s, Sadeh’s, Scripps’ and Oakley’s algorithms
Detection of consolidated periods of similar activity patterns
Crespo’s, Roenneberg’s algorithms
All algorithms have been implemented to return a binary time serie (0 being rest or activity depending on the definition made in the original article…)
In [1]:
import pyActigraphy
/usr/local/lib/python3.7/site-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.
import pandas.util.testing as tm
In [2]:
import plotly.graph_objs as go
In [3]:
import os
fpath = os.path.join(os.path.dirname(pyActigraphy.__file__),'tests/data/')
In [4]:
raw = pyActigraphy.io.read_raw_rpx(fpath+'test_sample.csv', period='6 days', language='FR')
All based on the same idea;
use a rolling window over the data and convolute them with a “gaussian”-like kernel.
if the resulting data is above a predefined threshold, classify as activity. Rest otherwise…
The weights of these “gaussian”-like kernels as well as the applied thresholds are very much device and population dependant. Using these algorithms would require to adapt these parameters to the population under study… Possible but seldom done.
In [5]:
layout = go.Layout(title="Rest/Activity detection",xaxis=dict(title="Date time"), yaxis=dict(title="Counts/period"), showlegend=False)
Disclaimer: in the original paper, they define “rescoring” rules… Not implemented in the package.
In [6]:
help(raw.CK)
Help on method CK in module pyActigraphy.sleep.scoring_base:
CK(scale=1e-05, window=array([ 400, 600, 300, 400, 1400, 500, 350, 0, 0], dtype=int32), threshold=1.0) method of pyActigraphy.io.rpx.rpx.RawRPX instance
Cole&Kripke algorithm for sleep-wake identification.
Algorithm for automatic sleep scoring based on wrist activity,
developped by Cole, Kripke et al [1]_.
Parameters
----------
scale: float, optional
Scale parameter P.
Default is 0.00001.
window: np.array
Array of weighting factors, :math:`W_{i}`.
Default is [400, 600, 300, 400, 1400, 500, 350, 0, 0]
threshold: float, optional
Threshold value for scoring sleep/wake.
Default is 1.0.
Returns
-------
ck: pandas.core.Series
Time series containing the `D` scores (0: sleep, 1: wake) for each
epoch.
Notes
-----
The output variable D of the CK algorithm is defined in [1]_ as:
.. math::
D = P*(
[W_{-4},\dots,W_{0},\dots,W_{+2}]
\cdot
[A_{-4},\dots,A_{0},\dots,A_{+2}])
with:
* D < 1 == sleep, D >= 1 == wake;
* P, scale factor;
* :math:`W_{0},W_{-1},W_{+1},\dots`, weighting factors for the present
minute, the previous minute, the following minute, etc.;
* :math:`A_{0},A_{-1},A_{+1},\dots`, activity scores for the present
minute, the previous minute, the following minute, etc.
References
----------
.. [1] Cole, R. J., Kripke, D. F., Gruen, W., Mullaney, D. J.,
& Gillin, J. C. (1992). Automatic Sleep/Wake Identification
From Wrist Activity. Sleep, 15(5), 461–469.
http://doi.org/10.1093/sleep/15.5.461
Examples
--------
In [7]:
CK = raw.CK()
In [8]:
layout.update(yaxis2=dict(title='Classification',overlaying='y',side='right'), showlegend=True);
In [9]:
go.Figure(data=[
go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
go.Scatter(x=CK.index.astype(str),y=CK, yaxis='y2', name='CK')
], layout=layout)
In [10]:
sadeh = raw.Sadeh()
scripps = raw.Scripps()
In [11]:
go.Figure(data=[
go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
go.Scatter(x=sadeh.index.astype(str),y=sadeh, yaxis='y2', name='Sadeh'),
go.Scatter(x=scripps.index.astype(str),y=scripps, yaxis='y2', name='Scripps')
], layout=layout)
This is the sleep/wake scoring algorithm used by the Actiware software from Respironics. The various activity thresholds have thus been tailored for Actiwatch devices. Be careful when applying this algorithm on data acquired with other devices.
In [12]:
help(raw.Oakley)
Help on method Oakley in module pyActigraphy.sleep.scoring_base:
Oakley(threshold=40) method of pyActigraphy.io.rpx.rpx.RawRPX instance
Oakley's algorithm for sleep/wake scoring.
Algorithm for automatic sleep/wake scoring based on wrist activity,
developed by Oakley [1]_.
Parameters
----------
threshold: float or str, optional
Threshold value for scoring sleep/wake. Can be set to "automatic"
(cf. Notes).
Default is 40.
Returns
-------
oakley: pandas.core.Series
Time series containing scores (1: sleep, 0: wake) for each
epoch.
Notes
-----
The output variable O of Oakley's algorithm is defined as:
.. math::
O = (
[W_{-2},W_{-1},W_{0},W_{+1},W_{+2}]
\cdot
[A_{-2},A_{-1},A_{0},A_{+1},A_{+2}])
with:
* O <= threshold == sleep, 0 > threshold == wake;
* :math:`W_{0},W_{-1},W_{+1},\dots`, weighting factors for the present
epoch, the previous epoch, the following epoch, etc.;
* :math:`A_{0},A_{-1},A_{+1},\dots`, activity scores for the present
epoch, the previous epoch, the following epoch, etc.
The current implementation of this algorithm follows the description
provided in the instruction manual of the Actiwatch Communication and
Sleep Analysis Software (Respironics, Inc.) [2]_ :
* 15-second sampling frequency:
.. math::
W_{-8} &= \ldots = W_{-5} = W_{+5} = \ldots = W_{+8} = 1/25 \\
W_{-4} &= \ldots = W_{-1} = W_{+1} = \ldots = W_{+4} = 1/5 \\
W_{0} &= 4
* 30-second sampling frequency:
.. math::
W_{-4} &= W_{-3} = W_{+3} = W_{+4} = 1/25 \\
W_{-2} &= W_{-1} = W_{+1} = W_{+2} = 1/5 \\
W_{0} &= 2
* 60-second sampling frequency:
.. math::
W_{-2} &= W_{+2} = 1/25 \\
W_{-1} &= W_{+1} = 1/5 \\
W_{0} &= 1
* 120-second sampling frequency:
.. math::
W_{-1} &= W_{+1} = 1/8 \\
W_{0} &= 1/2
The *Automatic Wake Threshold Value* calculation is this [2]_:
1. Sum the activity counts for all epochs of the data set.
2. Count the number of epochs scored as MOBILE for the data set
(the definition of MOBILE follows).
3. Compute the MOBILE TIME (number of epochs scored as MOBILE from
step 2 multiplied by the Epoch Length) in minutes.
4. Compute the Auto Wake Threshold = ((sum of activity counts from
step 1) divided by (MOBILE TIME from step 3)) multiplied by 0.88888.
*Definition of Mobile* [2]_:
An epoch is scored as MOBILE if the number of activity counts recorded
in that epoch is greater than or equal to the epoch length in 15-second
intervals. For example,there are four 15-second intervals for a
1-minute epoch length; hence, the activity value in an epoch must be
greater than, or equal to, four, to be scored as MOBILE.
References
----------
.. [1] Oakley, N.R. Validation with Polysomnography of the Sleepwatch
Sleep/Wake Scoring Algorithm Used by the Actiwatch Activity
Monitoring System; Technical Report; Mini-Mitter: Bend, OR, USA,
1997
.. [2] Instruction manual, Actiwatch Communication and Sleep Analysis
Software
(https://fccid.io/JIAAWR1/Users-Manual/USERS-MANUAL-1-920937)
A threshold of 40 corresponds to a medium wake threshold.
In [13]:
oakley = raw.Oakley(threshold=40)
It is also possible to compute a threshold value automatically, based on activity data.
In [14]:
oakley_auto = raw.Oakley(threshold='automatic')
In [15]:
go.Figure(data=[
go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
go.Scatter(x=oakley.index.astype(str),y=sadeh, yaxis='y2', name='Oakley (thr: medium)'),
go.Scatter(x=oakley_auto.index.astype(str),y=scripps, yaxis='y2', name='Oakley (thr: automatic)')
], layout=layout)
This is a threshold-based algorithm that used morphological filters to “clean” short periods of activity (rest) surrounded by periods of rest (acitivity).
In [16]:
help(raw.Crespo)
Help on method Crespo in module pyActigraphy.sleep.scoring_base:
Crespo(zeta=15, zeta_r=30, zeta_a=2, t=0.33, alpha='8h', beta='1h', estimate_zeta=False, seq_length_max=100, verbose=False) method of pyActigraphy.io.rpx.rpx.RawRPX instance
Crespo algorithm for activity/rest identification
Algorithm for automatic identification of activity-rest periods based
on actigraphy, developped by Crespo et al. [1]_.
Parameters
----------
zeta: int, optional
Maximum number of consecutive zeroes considered valid.
Default is 15.
zeta_r: int, optional
Maximum number of consecutive zeroes considered valid (rest).
Default is 30.
zeta_a: int, optional
Maximum number of consecutive zeroes considered valid (active).
Default is 2.
t: float, optional
Percentile for invalid zeroes.
Default is 0.33.
alpha: str, optional
Average hours of sleep per night.
Default is '8h'.
beta: str, optional
Length of the padding sequence used during the processing.
Default is '1h'.
estimate_zeta: bool, optional
If set to True, zeta values are estimated from the distribution of
ratios of the number of series of consecutive zeroes to
the number of series randomly chosen from the actigraphy data.
Default is False.
seq_length_max: int, optional
Maximal length of the aforementioned random series.
Default is 100.
verbose: bool, optional
If set to True, print the estimated values of zeta.
Default is False.
Returns
-------
crespo : pandas.core.Series
Time series containing the estimated periods of rest (0) and
activity (1).
References
----------
.. [1] Crespo, C., Aboy, M., Fernández, J. R., & Mojón, A. (2012).
Automatic identification of activity–rest periods based on
actigraphy. Medical & Biological Engineering & Computing, 50(4),
329–340. http://doi.org/10.1007/s11517-012-0875-y
Examples
--------
>>> import pyActigraphy
>>> rawAWD = pyActigraphy.io.read_raw_awd(fpath + 'SUBJECT_01.AWD')
>>> crespo = rawAWD.Crespo()
>>> crespo
2018-03-26 14:16:00 1
2018-03-26 14:17:00 0
2018-03-26 14:18:00 0
(...)
2018-04-06 08:22:00 0
2018-04-06 08:23:00 0
2018-04-06 08:24:00 1
Length: 15489, dtype: int64
In [17]:
crespo = raw.Crespo()
In [18]:
crespo_6h = raw.Crespo(alpha='6h')
In [19]:
crespo_zeta = raw.Crespo(estimate_zeta=True)
In [20]:
go.Figure(data=[
go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
go.Scatter(x=crespo.index.astype(str),y=crespo, yaxis='y2', name='Crespo'),
go.Scatter(x=crespo_6h.index.astype(str),y=crespo_6h, yaxis='y2', name='Crespo (6h)'),
go.Scatter(x=crespo_zeta.index.astype(str),y=crespo_zeta, yaxis='y2', name='Crespo (Automatic)')
], layout=layout)
Like for the other algorithms, the default parameters correspond to those described in the original papers, which, most likely, have been tuned to the population used to validate the algorithm.
Since the output is made of consolidated periods, it is possible to return the offset and onset times of each period:
In [21]:
aot = raw.Crespo_AoT()
In [22]:
aot
Out[22]:
(DatetimeIndex(['2015-02-04 12:45:00', '2015-02-05 06:13:00',
'2015-02-06 07:44:00', '2015-02-07 05:45:00',
'2015-02-08 06:22:00', '2015-02-09 07:09:00',
'2015-02-10 05:56:00', '2015-02-10 11:45:00'],
dtype='datetime64[ns]', freq=None),
DatetimeIndex(['2015-02-04 11:46:00', '2015-02-04 22:07:00',
'2015-02-05 19:49:00', '2015-02-06 19:43:00',
'2015-02-07 21:28:00', '2015-02-08 19:27:00',
'2015-02-09 19:42:00', '2015-02-10 08:44:00'],
dtype='datetime64[ns]', freq=None))
In [23]:
aot[0]-aot[1]
Out[23]:
TimedeltaIndex(['00:59:00', '08:06:00', '11:55:00', '10:02:00', '08:54:00',
'11:42:00', '10:14:00', '03:01:00'],
dtype='timedelta64[ns]', freq=None)
Also threshold-based but uses correlations with test series of various lengths to find the consolidated period that best matches best the data.
In [24]:
help(raw.Roenneberg)
Help on method Roenneberg in module pyActigraphy.sleep.scoring_base:
Roenneberg(trend_period='24h', min_trend_period='12h', threshold=0.15, min_seed_period='30Min', max_test_period='12h', n_succ=3) method of pyActigraphy.io.rpx.rpx.RawRPX instance
Automatic sleep detection.
Identification of consolidated sleep episodes using the
algorithm developped by Roenneberg et al. [1]_.
Parameters
----------
trend_period: str, optional
Time period of the rolling window used to extract the data trend.
Default is '24h'.
min_trend_period: str, optional
Minimum time period required for the rolling window to produce a
value. Values default to NaN otherwise.
Default is '12h'.
threshold: float, optional
Fraction of the trend to use as a threshold for sleep/wake
categorization.
Default is '0.15'
min_seed_period: str, optional
Minimum time period required to identify a potential sleep onset.
Default is '30Min'.
max_test_period : str, optional
Maximal period of the test series.
Default is '12h'
n_succ : int, optional
Number of successive elements to consider when searching for the
maximum correlation peak.
Returns
-------
rbg : pandas.core.Series
Time series containing the estimated periods of rest (1) and
activity (0).
References
----------
.. [1] Roenneberg, T., Keller, L. K., Fischer, D., Matera, J. L.,
Vetter, C., & Winnebeck, E. C. (2015). Human Activity and Rest
In Situ. In Methods in Enzymology (Vol. 552, pp. 257-283).
http://doi.org/10.1016/bs.mie.2014.11.028
Examples
--------
In [25]:
roenneberg = raw.Roenneberg()
roenneberg_thr = raw.Roenneberg(threshold=0.25, min_seed_period='15min')
/Users/ghammad/Work/CRC/pyActigraphy/pyActigraphy/sleep/scoring/utils.py:99: RuntimeWarning:
invalid value encountered in greater
In [26]:
go.Figure(data=[
go.Scatter(x=raw.data.index.astype(str),y=raw.data, name='Data'),
go.Scatter(x=roenneberg.index.astype(str),y=roenneberg, yaxis='y2', name='Roenneberg'),
go.Scatter(x=roenneberg_thr.index.astype(str),y=roenneberg_thr, yaxis='y2', name='Roenneberg (Thr:0.25)')
], layout=layout)
In [27]:
aot = raw.Roenneberg_AoT(threshold=0.25, min_seed_period='15min')
/Users/ghammad/Work/CRC/pyActigraphy/pyActigraphy/sleep/scoring/utils.py:99: RuntimeWarning:
invalid value encountered in greater
In [28]:
aot
Out[28]:
(DatetimeIndex(['2015-02-04 13:03:00', '2015-02-05 08:16:00',
'2015-02-05 10:56:00', '2015-02-05 15:08:00',
'2015-02-05 15:43:00', '2015-02-06 03:31:00',
'2015-02-06 07:36:00', '2015-02-06 11:15:00',
'2015-02-06 19:32:00', '2015-02-06 20:06:00',
'2015-02-07 07:47:00', '2015-02-07 11:01:00',
'2015-02-07 13:37:00', '2015-02-07 18:58:00',
'2015-02-08 07:35:00', '2015-02-08 14:31:00',
'2015-02-09 07:20:00', '2015-02-09 11:06:00',
'2015-02-10 06:39:00', '2015-02-10 07:32:00'],
dtype='datetime64[ns]', name='Date_Time', freq=None),
DatetimeIndex(['2015-02-04 12:45:00', '2015-02-04 21:07:00',
'2015-02-05 10:32:00', '2015-02-05 14:41:00',
'2015-02-05 15:24:00', '2015-02-05 19:25:00',
'2015-02-06 03:50:00', '2015-02-06 10:56:00',
'2015-02-06 18:58:00', '2015-02-06 19:44:00',
'2015-02-06 20:37:00', '2015-02-07 10:42:00',
'2015-02-07 13:19:00', '2015-02-07 18:36:00',
'2015-02-07 21:00:00', '2015-02-08 14:01:00',
'2015-02-08 20:10:00', '2015-02-09 10:42:00',
'2015-02-09 18:55:00', '2015-02-10 06:48:00'],
dtype='datetime64[ns]', name='Date_Time', freq=None))
In [29]:
aot[0]-aot[1]
Out[29]:
TimedeltaIndex(['00:18:00', '11:09:00', '00:24:00', '00:27:00', '00:19:00',
'08:06:00', '03:46:00', '00:19:00', '00:34:00', '00:22:00',
'11:10:00', '00:19:00', '00:18:00', '00:22:00', '10:35:00',
'00:30:00', '11:10:00', '00:24:00', '11:44:00', '00:44:00'],
dtype='timedelta64[ns]', name='Date_Time', freq=None)
Et voilà! Easy, isn’t it?